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Abstract 



The critical behaviour of a bond-disordered Ashkin- Teher model on a 
square lattice is investigated by intensive Monte-Carlo simulations. A duality 
transformation is used to locate a critical plane of the disordered model. This 
critical plane corresponds to the line of critical points of the pure model, 



Tlj" I along which critical exponents vary continuously. Along this line the scaling 

^-^ ' exponent corresponding to randomness (p = (ajv') varies continuously and 

is positive so that randomness is relevant and different critical behaviour 
^ ■ is expected for the disordered model. We use a cluster algorithm for the 

Monte Carlo simulations based on the Wolff embedding idea, and perform a 
C^ \ finite size scaling study of several critical models, extrapolating between the 

critical bond-disordered Ising and bond-disordered four state Potts models. 
'^ ■ The critical behaviour of the disordered model is compared with the critical 

^ ■ behaviour of an anisotropic Ashkin- Teller model which is used as a refference 

Q . pure model. We find no essential change in the order parameters' critical 

exponents with respect to those of the pure model. The divergence of the 
k>( ^ specific heat C is changed dramatically. Our results favor a logarithmic type 

Vh ' divergence at T^^ C ~ log L for the random bond Ashkin- Teller and four state 

Potts models and C ~ log log L for the random bond Ising model. 
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I. INTRODUCTION 

How is the critical behaviour affected by the introduction of disorder ( usualy dilution or 
bond-randomness ) into a model ? The Harris criterion [|l]-^ states that (j), the scaling index 
of the operator corresponding to randomness at the pure system fixed point (also called 
the crossover exponent ) is equal to - {a and u are the specific-heat and correlation length 
exponents of the pure model). Thus the critical behaviour of the pure system is unaltered 
by disorder if a < 0, and altered when a > 0. If a = the situation is marginal. 

Renormalization-group methods (namely expansion in the parameter e = 4 — d, where d 
is the dimensionality), applied to n-component continuous spin models with weak quenched 
bond disorder 0], confirmed the Harris criterion (p = -. Moreover, it was found that for 
n < ric = 4 — 4e + 0(e^) a new stable 'random' fixed point with new exponents appears. For 
the case n = 1 (Ising model) the stability of an 0(e^/^) random fixed point was shown by 
Jayaprakash and Katz [Q. For the d = 2 Ising model a = 0, and the operator corresponding 
to randomness is expected to be marginally relevant. 

Two partially conflicting theories by Dotsenko and Dotsenko (DD) and by Shalaev 
0, Shankar 0, and Ludwig ^ (SSL) were suggested for the d = 2 Ising model with small 
randomness. While both theories agreed on the divergence of the specific heat, e.g. 

Cimp(t)~lnlnl/|t| , (1) 

conflicting predictions were made for other quantities like the magnetization M and sus- 
ceptibility X- Extensive Monte-Carlo (MC) simulations of dilute and random-bond Ising 
models ( see a review in ref. [^ and references therein) helped to confirm the prediction 
on the specific heat. The other quantities were seen to behave ( as predicted by SSL), in 
essence, as in the pure model. The disordered d = 3 Ising model has been studied by means 
of MC simulations as well |T^ . 

Less attention has been paid in recent times to other disordered systems such as random 
q-state Potts models. Kinzel and Domany used a real-space Migdal-Kadanoff renor- 
malization transformation for the random-bond Potts model on the square lattice in two 
dimensions and found ( when apure > 0) a stable random fixed point with a reduced (with 
respect to the pure model) but positive arandom- Using methods of the same spirit, Andel- 



man and Berker |jTT[] also found a random fixed point. For g = 4 their calculation suggests 



a = —.37 and u = 1.19. Other workers ||r^ , p!3[| too have suggested a negative a for q > 2 



Potts models in two dimensions. MC work on random q-state Potts models has not been car- 
ried out, except for recent work on the bond-random two dimensional 8-state Potts model 



I^ , where randomness was shown to change the first order phase transition of the pure 



model into a second order one [|T5|. Novotny and Landau (NL) [0 used MC simulations to 



study the two dimensional site-dilute Baxter- Wu |]T^ model which is of the same universal- 



ity class as the 4-state Potts model. NL measured the exponents a = 0.0(2), u = 1.00(7), 
and 7 = 1.95(8) which did not seem to depend on the amount of dilution. They could not 
determine whether the new stable fixed point is a d = 2 random Ising fixed point (which is 
very similar to the pure Ising fixed point) or a new impure Baxter- Wu fixed point. 

In this paper we address the effect of randomness on a critical line which connects the 
critical points of the d = 2 Ising and 4-state Potts models. Such is the critical line of the 
two dimensional Ashkin- Teller model |TB[ on a square lattice. 



A convenient representation of the pure Ashkin- Teller (AT) model is in terms of two 
Ising spin variables, ai and Xj, placed on every site of a lattice. Denoting by < ij > a pair 
of nearest neighbor sites, the Hamiltonian is given by 

n = - Y^ [KaiCj + KTiTj + AaiajTiTj] (2) 

Here K is the strength of the interactions between neighboring a and r spins, and A is a 
four-spin coupling ( throughout this article we absorb a factor of l/ksT into the coupling 
constants). The phase diagram of the ferromagnetic AT model is known in two dimensions 
from duality transformations and renormalization group studies [0,^. The three dimen- 



sional model has been studied as well |^ . The phase diagram which we reviewed previously 
2^ includes a line of critical points which connects the Ising critical point at one end to 
the 4-state Potts critical point at the other end. Along this line the exponents vary continu- 
ously, and have been determined analytically p3|j2^ , interpolating between their Ising and 



four-state Potts values. For instance, the crossover exponent connected with randomness 
= ^ changes smoothly from at the decoupled Ising point, to 1 at the 4-state Potts point. 
Therefore, according to the Harris criterion, we expect randomness to be a relevant operator 
of varying strength in this regime of the AT model, and expect the critical behaviour of the 
disordered model to differ from that of the pure model. 

The critical behaviour of the disordered AT model has not been studied before, apart 



from a conjecture by Alcaraz and Tsallis [^ as to the location of the critical manifold of the 



bond dilute AT model. None the less the effect of disorder on the critical line of the Baxter( 



or symmetric 8- vertex) model p6|], which is isomorphic to the critical line of the AT model 
p3| , |2^ , was considered. DD [^ have extended their study of the disordered Ising model to 
the Baxter model, and found, for small disorder, the specific heat divergence to change from 
its pure form to the form of eq. (0). Matthews- Morgan, Landau and Swendsen (MLS) ||2^ 
have studied the site-dilute Baxter model by means of a MCRG ( MC renormalization-group) 
method. MLS found - = yr = 0.98(7) and - slightly below the Ising value. Both studies 
indicate that continuous variation of critical exponents of the pure model is substituted by 
flow to a random Ising fixed point. 



In a previous article |2^ we described an efficient MC cluster algorithm |2^ for the 
AT model. In this work we have used this algorithm to perform extensive MC simulations 
of the random-bond AT model. The importance of using a cluster algorithm for the 



acquirement of reasonably accurate data cannot be overlooked. This is especially true in the 



highly random regime where local algorithms are extremely slow |^ while cluster algorithms 
are most efficient [p2| , p3| .The partial elimination of critical slowing down allows us to go to 
large lattice sizes (up to L=256) and obtain accurate results on the finite size scaling of 
the thermodynamic quantities at criticality. Our results seem to favor an Ising like critical 
behaviour. 

This work is organized as follows. In sec. ^ we define the random-bond AT model 
and a related anisotropic AT model. As will be explained, the anisotropic AT model is 
used as a non-disordered reference model for the random-bond model. We use the duality 
transformation of the AT model to locate exactly critical surfaces of the two models, which 
are related to the line of critical points of the pure model. In sec. |T| we describe the 
MC cluster algorithm, give definitions of the measured quantities and describe our choice of 



simulation points. We have performed simulations of the anisotropic and random-bond AT 



models at several points on their critical surface. In sec. |I^ we give our results which consist 
of the finite size dependence of thermodynamic quantities of the two models at criticality. 
Analysis of our results shows that the anisotropic model exhibits varying critical exponents, 
i.e. the specific heat diverges as 

C(t)~|t|-" , (3) 

where at different points of the critical surface of the anisotropic AT model we measured 
different values of a. On the other hand, the random-bond model seems to exhibit a single 
type of behaviour, i.e. the specific heat diverges as 

amp(t)~-ln|t| . (4) 

this is identical with a pure Ising behaviour and not with a random Ising one. In sec. we 
summarize our results. 



II. DEFINITION OF MODELS AND LOCATION OF CRITICAL PLANES 

The model we wish to study is the Random-Bond AT model (RBAT) on a square lattice, 
which is defined by the Hamiltonian 

7Y = - ^ [Kijaiaj + KijTiTj + AijaiT.ajTj] . (5) 

The coupling constants Kij and Aj ,,• are chosen according to 

{z<^ A \ — i (-^^'^^) with probability p , . 

lA^j, A,jj = j (^2^ ^2) ^i^i^ probability 1-p ' ^^^ 

It seems quite natural ( the need for this comparison will be made clear at the end of 
this section) to compare the RBAT with the Anisotropic AT model (AAT) with the same 
Hamiltonian (^ but with the couplings distributed as follows: 

(T^ A \ — j (-^^^^) fo^ bonds (z,j) in the horizontal direction ,„. 

^ *'■'' ^J'^ ~ j (i^2,A^) for bonds (i,j) in the vertical direction ' ^ ^ 

Spatial anisotropy is usualy known to be an irrelevant operator, that is, the critical 
behaviour of a model is unaltered by introduction of anisotropy. Thus one expects that 
continuous variation of critical exponents will continue to exist on some critical manifold of 
the AAT, as is the case for the isotropic (or pure) model. 

A duality transformation can be used to locate critical manifolds of the AAT and RBAT 
in a way we now present. 



A. Duality and location of a phase transition 

Consider an AT system on a square lattice with spatially varying coupling constants (the 
spatial variation can be due to anisotropy, or it can be due to randomness as we consider 
later). The dual lattice D of the original square lattice G is the square lattice whose sites are 
at the centers of the plaquettes of the original lattice. In figure |l| sites of G are denoted by 
dots, while sites of D are denoted by crosses. Consider a bond of strength {K, A) between 
two sites in G. Under duality it is mapped onto a bond of strength {K^K) = D[{K,A)], 



where D is the duality transformation of the AT model ||3^ , |20[| which we list in the Appendix. 
The dual bond lies in the dual lattice D, intersecting the original bond. This situation is 
depicted in figure |l|. The original bond of strength {K, A) is depicted by the solid line while 
its dual is depicted by the dashed line. 

Suppose the only spatial variation of the couplings is one of anisotropy. Then an AT 
system with couplings of strength {K^,A^) in the horizontal direction and strength {K"^, A^) 
in the vertical direction, transforms under duality into an AT system on the dual lattice D 
with couplings of strength {K^,A^) in the vertical direction and strength {K'^,A'^) in the 
horizontal direction. The free energies of the original system and its dual have the same 
singular part, and therefore we can write 

U{K\ A'), {K\ A')] = a{K\ A^), {K\ A^)] , (8) 

where the first pair of couplings corresponds to the horizontal bonds, and the second pair 
to the vertical ones. The duality transformation D is self-inverse. Thus,jin anisotropic 
AT system with couplings [(/f^ A^), (i^^, A^)] such that {K^^A"^) = {K'^,A'^), transforms 
under duality onto itself. So, there exists a two dimensional self-dual manifold (it would 
help the reader to view the duality transformation as a mathematical mapping in the four 
dimensional space of couplings [{K^,A^), (K^, A^)] ). 

The line of critical points of the isotropic AT model lies within the subspace S, defined 
by < A < Ji'. E is located between a decoupled Ising model where A = and a four-state 
Potts model with A = K. Two sub-areas of S which map onto each other under duality are 
separated by a self dual line ( fulfilling {K,A) = {K,A) ). The assumption that S contains 
only two phases, with a single phase transition between them(where fs is non-analytic) leads 
to the conclusion that this phase transition must occur on the self-dual line. In this way the 
line of critical points of the isotropic AT was located exactly [pDj]. 

Let us now consider the corresponding subspace, S^, of the AAT for which < A^ < 
K^ and < A^ < K^. Ea lies between an anisotropic decoupled Ising subspace where 
A^ = A^ = and an anisotropic four-state Potts subspace with A^ = K^, A^ = K^. S^ is 
invariant under duality since for Q < A < K we get (see Appendix) Q < A< K. It follows 
also that the self dual plane < A^ < K^ , (K^, A^) = {K^^A^) is included in the subspace 
we are considering. Thus, again assuming that S^ contains only two phases, with a single 
phase transition between them, we conjecture that the self dual plane is a plane of phase 
transitions (this kind of argument was first used by Fisch |^ for the random-bond Ising 
model, and later for Potts models by Kinzel and Domany |^). Since the space of couplings 
is four dimensional, while the self dual plane is only two dimensional, the self dual plane is 
only part of the three dimensional critical manifold in S^ . Notice that this is unlike the 



case of the Potts model or any model with a single coupling constant where the self dual 
line is identical with the critical line. 



B. A Critical Plane of the RBAT 

The argument that leads to locating the corresponding critical plane of the RBAT is 
very similar to the one presented above, up to a subtle difference. Consider a specific 
realization of a RBAT on a square lattice with couplings {K^,A^), (K^, A^) distributed at 
random with probabilities p and 1 —p, respectively. Replacing each bond by its dual we get 
(on the dual lattice) a system with couplings {K^,A^) and (i^^, A^), distributed at random 
with probabilities p and 1 — p. This defines a one to one mapping under duality from the 
ensemble of realizations of [{K^,A^), {K'^,A'^),p] RBAT systems, onto the ensemble of the 
dual systems [{K^,A^), {K'^,A'^),p]. Under duality each bond is mapped onto the bond on 
the dual lattice which intersects the original bond (see fig. |l|). Thus, the spatial distribution 
of the two types of bonds on the dual lattice is different from that of the original lattice. 
But since all spatial distributions have the same weight in their ensemble, there is a one to 
one correspondence of bond configurations in the two ensembles, and we can write: 

U{K\A'),{K\A'),p] = a{K\M),{IO,A^),p], (9) 

where it is understood that in calculating /^ we average over all spatial bond distributions. 
Choosing p = | we find that there is a self dual manifold at {K^, A^) = {K^, A^). Our ability 
to locate a self dual manifold for the case p = | makes this choice of p a convenient case for 
studying the RBAT. Viewing the duality transformation as a mathematical mapping in the 
space of couplings [{K^, A^), {K^, A^), |], it is identical to the duality transformation in the 
space of couplings [{K'^,A^), (i^^,A^)] of the anisotropic AT model. Thus all conclusions 
which were based on duality in sec. |11A| apply here equally. Following the logic of the 



previous subsection we conjecture that the subspace of the RBAT for which < A^ < K^ 
, < A^ < K^ and p = \ includes in it a critical plane defined by < A^ < K^ , 
{K\A^) = (k\A^) . 

In order to study the effect of randomness on the AT model we need a reference non- 
disordered system for comparison. When studying a specific RBAT model with couplings 
[{K^,A^), {K^,A'^),^] it is not clear with which specific non-random AT model with cou- 
plings {K, A) should one compare. None the less a comparison with the corresponding 
AAT model with couplings [(i^^,A^), (i^^,A^)] seems to be quite natural. Consider the 
anisotropic model with bonds {K^,A^) in the horizontal direction and bonds (i^^,A^) in 
the vertical direction. Next, consider as a perturbation, changing a small fraction p of the 
bonds in the horizontal direction into bonds of strength {K'^,A'^) and the same fraction of 
bonds in the vertical direction into bonds of strength (if^,A^). When p increases all the 
way to p = ^ the anisotropy disappears and the model is completely random. So, in a 
sense, the perturbation of randomness can take us continuously from the critical plane of 
the anisotropic model to the critical plane of the random model. In view of this simple 
picture we have performed simulations of the anisotropic and random-bond AT models at 
several points on their self-dual planes described above. We indeed find the two planes to 
be critical. In the next section we give some details of the simulations. 



III. METHOD AND DETAILS OF SIMULATIONS 

A. The Method 

In ref. [^ we described a cluster algorithm for the AT model and tested its performance 
(on the isotropic model), as explained there, the algorithm can be regarded from two 
different points of view. The conceptually simpler and faster ( though less general ) one is 
that of the Wolff embedding idea |]36| , |37| . It is faster in the sense that its implementation 
on the computer runs faster since the freeze-delete procedure is less complicated. The main 
idea is to embed into the AT model an Ising model of space dependent couplings Jjk and 
simulate it using the SW |^ or Wolff |^ procedure for the Ising model. To be explicit, 



consider the Hamiltonian (^, and take the r variables as fixed, so we can write 

U = Hi+H2 = - Yl (^iJ + K,jnrj)criaj - ^ KijTiTj. (10) 

7^2 represented by the second sum is a constant, and remembering that we are considering 
Ajj < Kij for all (i,j), Hi is a ferromagnetic Ising model in the a variables with couplings 
Jij = Kij + AijTiTj. Simulating this Ising model with any procedure that will maintain 
detailed balance with respect to TYi and will not change the value of 7^2, will also maintain 
detailed balance with respect to Ti. So we can use for example the SW pSf or Wolff iSH 



procedure for the Ising model. This by itself will maintain detailed balance but will not 
be ergodic since the r variables will not be updated. Obviously to update the r variables 
the same process should be repeated , holding the a variables fixed and simulating an Ising 
Hamiltonian with the r variables. To summarize, our procedure goes as follows: Choose at 
random whether to embed into the AT Hamiltonian an Ising Hamiltonian in the a (or r) 
spins. Pick a random site in the lattice, grow a cluster of a (or r) spins, using the Wolff 
single cluster procedure with the Ising Hamiltonian ( ( PUJ ) or the opposite one in case of a r 
embedding), and flip it. 

B. Choice of Simulation Points 

All our simulations were performed at the critical plane defined by 

{[{K\A'),{K\A'),p=^]: {K\A') = {K\A^) , < A' < K' } , (11) 

as explained in sec. H (for the anisotropic model just omit the p = |). For convenience we 
also use the following notation for the coupling constants: 

Z = exp-4^ ; X = exp-2(^^'+^) ; X = {X,Z) . (12) 

The first series of measurements were performed at five points {X^i,X'^i), i = 0. . .4, 
which we label as Cj, i = . . A. These were chosen so as to interpolate between {X^q, X'^q), 
which is an anisotropic (or random-bond) decoupled Ising critical point, and {X^4,X'^4), 
which is an anisotropic (or random-bond) four-state Potts critical point. The points Cj 
interpolate in a similar manner to the way in which the critical line of the pure AT connects 



the pure decoupled Ising critical point with the pure four-state Potts critical point. Thus 
we chose the five points X^i to lie equidistantly on the fine 

Z - ^Ising = -2(X - Xjgj^g) . (13) 

The slope —2 was chosen so that the line ([l3|) lies parallel (in the XZ plane) to the critical 
line of the pure AT model, Z = 1 — 2X. Note that once {X^, Z^) are defined, (X^, Z^) are 
defined through (|ll], p^ . 

We clearly have some freedom in our choice of (Xj^sing^ ^Lng); the decoupled Ising cou- 
plings; setting A^ = in (|ll|) one gets A^ = A^ = as well. The remaining parameter, K^, 
is chosen so that ( see (pjT) ) 
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Ising Ising in Ising" \ / 



Since the extent of deviation from pure behavior is obviously determined not only by p but 
also by the difference between the two sets of couplings, the ratio of ^ was chosen so that 
randomness will be pronounced |]32|j39| . 



While (14) defines the first simulation point X^q = X^isingj the fifth simulation point 
X^4 = X^Potts is chosen to be a four-state Potts point lying at the intersection of the line 
(H) with the four-state Potts line Z^ = X^ (or A^ = K^). 

The other three simulation points X^i 2,3 lie equidistantly on the line (|T3|) between X^ising 
and XI Potts : 

^i = i(^Potts - ^liing) + ^liing ^i = ^Ling " 2(X/ - X^-^^) i = . . . 4 , (15) 

Thus the first series of points Ci at which we made measurements Co.. .4 is defined through 
equations (|11],|T2|,|T^|T5|) and the definition of X^Potts- 

Two additional measurement points were intended to monitor the effect of varying the 
amount of randomness in the model on the critical behaviour. In order to define these two 
points we define two series of five points each, A0...4 and -B0...4 • Similarly to the series Co.. .4, 
these two series are defined through equations (|Tl|JI^Jl5|) . But in contrast with the series 
Co.. .4, eq. (|T3D is substituted by 



for the series A0...4, and by 



7^2 1 

Ising -■- 

-falsing ^ 



K'^ 1 

Ising -■- 

Ising 



(16) 



(17) 



for the series -Bo.. .4- Though we defined here ten additional points we actually performed 
measurements only at the points A2 and B2. The points A2, B2, C2 represent three RBAT 
( or AAT ) models with coupling ratios ^ ~ | and -^ ~ |, \i jq respectively. 



C. MC Simulations 



We calculated the energy [^ per site: 

1 



E = [{E)] 



L2 



[( Y. {^iA^i^J + ^i^j) + Kj'^i^i'^jTj})] 



(18) 



<»,j> 



where L is the linear lattice size; the angular brackets denote the usual thermal MC average 
whereas the square brackets denote the quenched average over random-bond configurations 
. The specific heat per site follows from the energy fiuctuations 



C 
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We define 

To compute the magnetization per site, M, we used 



M = -[{\m^\ + \m^\)] 



12 J2i O-iT 



(19) 



(20) 



(21) 



taking advantage of the symmetry between the a and r spins to increase accuracy. Another 
order parameter of the transition plane we are considering is the polarization [^], P, defined 

as 

^=[(K.|)], (22) 

which measures the magnetization of the err spins. We calculated the magnetic susceptibility 



of the a spins using [^ 



X= -ir[{ml + ml)] 



(23) 



and also measured c, the size of the cluster fiipped at each step, and calculated (c). There is 
a connection between the size of the clusters and the susceptibility |^ which is common for 
all algorithms that generate non-interacting clusters of spins (if all spins in a cluster have 
the same value). For the Wolff single cluster procedure it has the simple form [Q 



X= (c). 
Lastly we measured the susceptibility of the ar spins, Xp'- 

Xp = L'[{ml)] 



(24) 



(25) 



For each sample we estimated the effective autocorrelation times [^ of the energy, 
r|f and the magnetization, rff using the method outlined in ref. |3^. The length of an 
individual run was often too short for a precise estimate, but averaging over all samples 
gave a rough estimate of [t'^m] ■ These estimates of r^^ allowed us to verify that all samples 



9 



were equilibrated well enough before we started collecting data. Typically we discarded at 
least 50r|f spin configurations. 

The fiuctuations in the random bond configurations introduce, in addition to the usual 
MC error, another source of statistical error. Thus the total error for the specific heat, 5C, 
for example, is given by 



.2 ^' ^' 



^SCy = ^ + ^' , u.Tmc large. (26) 

where n is the number of random bond samples, Tmc is the length of the MC runs and r is 
the autocorrelation time of the algorithm ( which we can estimate by t^^). a^ is the variance 
of C within a given sample and is due to thermal fiuctuations, while a,^ is the variance of the 
exact (C) within the ensemble of random bond samples. Thus both the thermal fiuctuations 
and the sample to sample fiuctuations contribute to the error. 

We found that both terms in ( ]26|) imply that simulations of the AT models and the Potts 
model are much more expensive in computer time than simulations of the Ising model. First 



there is an increase in r in consistence with our findings in the pure models |2^. E.g. at 
L = 256 we find for the random Ising model t^ ~ 3-1(1) compared with t^ ?» 11.3(3) for 
the random four-state Potts model. Second we found that for all thermodynamic quantities 
P the relative variance cr^/-P^ was larger for the random four-state Potts model than for 
the random Ising model. This effect was strongest for the specific heat, e.g. for L = 256 
a^ c/C^ of the random four-state Potts model was approximately an order of magnitude 
larger than that of the random Ising model. In accordance with these findings, the number 
of samples simulated for the Ising model varied between n = 1000 for L = 4 and n = 150 
for L = 256, while for the four-state Potts model it varied between n = 2000 for L = 4 and 
n = 330 for L = 256. 



IV. MONTE CARLO RESULTS 

A. Specific Heat 

We expect the critical specific heat of the AAT to follow the finite size scaling [Q form 
of(i, 

C^ai + 61L"/'' . (27) 



The exponent ratio - is predicted analytically ||23|j2^ to vary continuously along the critical 
line of the pure AT model between - = at the decoupled Ising critical point and - = 
1 at the four-state Potts critical point. One expects the same variation to hold for the 
corresponding critical manifold of the AAT at which the simulations were performed. In 
figure H we plot the critical specific heat, C as a function of logL for the five critical AAT 
models Co. a- Fits are made to the form (^), and the estimated values of a/u are listed in 
Table |I[ Since the model Co is an anisotropic Ising model, its specific heat does not follow 
the form (P^), but rather C ~ A + B log L. The anisotropic four-state Potts model C4 should 
have a/u = 1 with a logarithmic correction [|^ to scaling C ^ L/ log^'^ L . In the range 



of L considered, the logarithmic correction results in an estimate for a/p which is too low. 

10 



The results for the specific heat of the AAT models clearly demonstrate the variation of ajv 
from its Ising value, ajv = to its four state Potts value. As will be seen shortly this is 
quite different from the results for the RBAT models. 

In order to analyze the specific heat results of the RBAT we first tried to fit the results 
to (0), but unlike the specific heat of the AAT, this was not possible for the full lattice 
size range 4 < L < 256. Thus in fig. |^ we show the specific heat as a function of log L for 
the five critical RBAT models Co. .4 along with fitting curves to (|27|), where only data for 
L > 24 were used to fit. The estimated exponents ajv obtained from these fits are within 
errors indicating a possible logarithmic behaviour. In addition, fig. ^ indicates a crossover 
of the effective exponents (the slope between every two successive data points) from higher 
to lower values. Physically we expect that this is due to a crossover from the pure models' 
exponents to some new random critical behaviour. We are thus led to the following finite 
size scaling ansatz for the specific heat 

C = ao + 60 ln[l + Co(L(°/'^)p - 1)] , (28) 

where [ajv^ is the critical exponent ratio of the corresponding pure model and cq = . ?■, . 
An important point to note is that for each random bond AT model, following the discussion 
at the end of sec. |IIB| , (a/z/)p is taken from our results (listed in Table |[) for the corre- 
sponding anisotropic AT model. (? is a measure of the randomness, and is related through 
the relation 

Co = (L("/^)- - l)-i , (29) 

to a crossover length L^, at which crossover from the pure model's power law behaviour to 
the random logarithmic behaviour occurs. Thus for L <^ L^. our eq. ( P5| ) reduces to ( P?] ) 
while for L ^ L,, and L^°'/'^^p ^ 1 a logarithmic behaviour is attained, 

C = a + b\nL. (30) 



Apart from crossing over to the correct pure result (|27|) when cq — *> 0, eq. (|28 
constructed so that in the Ising model limit, {a/h')p — > 0, it becomes 



was 



C = a + 61n(l + ^lnL) . (31) 



This is the finite size scaling form of (0) which was shown by Wang et al. |^2[ to fit the 
random bond Ising model well. This form will be examined below as another candidate for 
describing the RBAT models. 

In fig. ^ the specific heat of the five critical RBAT models C0..4 are plotted again but 
with fits to (^) using the full lattice size range 4 < L < 256, on a semilogarithmic scale. 
The fitting parameters are listed in Table |l| . For these points ( with large randomness, 
%T ^ jq ), the crossover lengths L^ are found to be 1. 

In fig. § the specific heats of three critical RBAT models A2,B2, C2 are plotted with the 
fitting functions of (p8|). (a/z/)p is taken from the corresponding anisotropic models and the 
fitting parameters are listed in Table H . (a/z/)p of the three models ^2,-82, C2 is very similar 
( .40, .37, .37 respectively) but they differ in their amount of randomness, %t ^ ^, Z' jq 
respectively. The form ( pSD seems to describe all three models adequately with the main 
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change being in the crossover length L^ which decreases as randomness increases. The 



estimates for Lc can be compared with results from random bond Ising models |]T0|,^ with 
the same coupling ratios. For ^ ~ ^, \) ^ we obtain Lc = 1, Lc = 4. ± .4, Lc = 51 ± 7 
compared with Lc = 2 ± 1, Lc = 16 ± 4 and L^. ~ 1000 for the respective Ising models. 
According to the Harris criterion one expects the crossover length to scale roughly as 

Lc ~ R'-^ , (32) 

where = {a/iy)p is the crossover exponent and i? (i? < 1, i? ~ ^)is a measure of 
the randomness. For the RBAT models considered here (a/z/)p ^ .37, while for the random 
bond Ising models randomness is marginal ( (a/i/)p = 0), so that the much smaller crossover 
lengths that we find are consistent with the Harris criterion. 

Since there is no theoretical prediction for the critical behaviour of the RBAT model, we 
examined possible scaling forms for the specific heat, other than (|28|) . A natural candidate 
is a double logarithmic behaviour as in eq. (|l|). This form was predicted by DD for the 



disordered Baxter |2^ and Ising |1T0|J3^ models and confirmed for the latter model. In its 



asymptotic limit, the finite size scaling form of (|I]) is 

C = a + 61nlnL. (33) 

We compare this form with (RO) as both these forms are asymptotic and do not include 



crossover from the pure behaviour. In figure ^ we plot again our specific heat results for all 
the models tested except for A2 which is not considered due to its large crossover length. 
The results are fitted by ( pOD in fig. || and by (^) in fig. |^, for L > 16. 

As found by Wang et al., for the random bond Ising model Cq, agreement with ( p3D is 
much better than with (|30D, according to our error analysis. For the model Ci too agreement 
with (|33D is better but this can be attributed to slow crossover behaviour from the small 
power law (a/z/)p ?^ .17 ( the condition L^^'/'^^p ^ 1 is not fulfilled for such small {a/iy)p and 
for the range of L considered). For the models 6*2,3 the quality of both fits is comparable 
while for the random bond four state Potts C4 the fit (|30|) is better. For the model B2 we 
have obtained larger statistics ( more than 900 samples for L < 128 and 300 samples for 
L = 256) and the result is conclusively in favor of (^). 

Last we considered the possibility that the crossover from the power law behaviour to a 
logarithmic behaviour is yet followed at even larger L by a crossover to a double logarithmic 
behaviour. This was done by fitting the specific heat data for L > 16 according to (pij). 
Though the fits were of good quality, the estimated crossover lengths Lj = exp{l/g) did 
not make sense physically. For the B2 model we found Li > 10^ which actually implies a 
logarithmic behaviour for all accessible values of L. This should be compared against our 
estimate Lc ~ 4 and the value Lj ^ 16 for the corresponding Ising model (see above). For 
the models Co.. .4 we find that Lj increases monotonically from the Ising model Cq where 
Lj Ri 2 to the four state Potts model C4 where Lj ?» 200. This is contrary to the expectation 
that the crossover length should decrease due to the increase in the crossover exponent 



(a/z/)p, as implied from (|3^) 



We conclude that our results are in favor of a logarithmic divergence of the specific 
heat and not a double logarithmic one for the RBAT model, excluding the random bond 
Ising model. The finite size scaling form ( pSD gives the best agreement with our results for 
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all the models, including those with weak randomness, and for the full lattice size range 
4 < L < 256. Thus we predict for C{t) the scaling form: 

C ~ ln[l + 6(t'"p - 1)] , (34) 



A possibility which we tend to rule out is that both types of critical behaviors ( that is ( |3^ 
and (|l]) ) occur and that there is some bifurcation plane which separates the two types of 
behaviors. 



B. Susceptibility and Magnetization 

1. Susceptibility 
According to finite size scaling theory one expects that the critical susceptibility diverges 



as 



X^ L- (35) 

for large enough L. The exponents' ratio 7/z/ is predicted analytically p3^p4| to be constant, 
■y/u = |, all along the critical line of the pure AT model. This is expected also to hold 
for the corresponding critical manifold of the AAT. Figure ^ shows results for the critical 
susceptibility as a function of logL for the five critical models, Co. a, of the AAT. So that 
the points won't fall on top of each other, x for the model Ci has been multiplied by 2*. The 
solid lines are linear fits to the form ( pSj ) for L > 24, yielding estimates for critical exponents 
^ which are listed in Table |. All results seem to be consistent with the analytic prediction, 
giving a combined best estimate of 7/z/ = 1.7512(6). An effective exponent analysis (or a 
trend analysis) shows that as logL increases, (^)eflF approaches the value | from above. 

The finite size scaling form ( |5BD and the value of ^ = | have been predicted []^-|^ for 
the disordered decoupled Ising model and confirmed in MC simulations [|^,0. Figure |^ 
shows results for the critical susceptibility as a function of log L for the five random critical 
models Co.. 4, of the RBAT. The same analysis has been carried out as for the AAT and the 



estimated critical exponents - are listed in Table |lTl[ The estimates for the decoupled Ising 
critical point Cq and for the models Ci^2 are consistent with ^ = |. On the other hand for 
the model C3 and the random-bond four-state Potts model C4 we obtain ^ = 1.736(3) and 
^ = 1.714(5) respectively. We also tried fitting the results with the form 

X = A^L^\nL)-P-, (36) 

for L > 24. The obtained estimate for py. is listed in Table |T|. for the random-bond four 
state Potts we find p^ = .148(21) while for the models C0-2 the error in p^ is at least of the 
same order as the estimate for p^ itself. 

2. Magnetization 



The critical magnetization of the five anisotropic models is expected [0,0 to follow the 
form 
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M ^ L-^ , (37) 

with f = |- Applying straightforward finite size scahng analysis, as described above for the 
susceptibility, to the magnetization results, we obtained estimates for - which are listed in 
Table U The results are slightly below the expected value, but an effective exponent analysis 
(or a trend analysis) shows that as logL increases (^)eff approaches the value | from below. 
The finite size scaling form (|37D and the value of ^ = | have been predicted [0-§| for 
the disordered decoupled Ising model and confirmed in MC simulations jT^JS^. Figure PUJ 
shows results for the critical magnetization as a function of log L for the five critical models, 
Co. .4, of the RBAT. A fit to the form (pTf ) yields the estimated critical exponents - listed in 



Table |T|. The estimates for the decoupled Ising critical point Cq and for the models Ci^2 are 



consistent with - = k- On the other hand for the model C3 and the random-bond four-state 
Potts model C4 we obtain - = .133(2) and - = .145(3) respectively. We also tried fitting 
the results with the form 

M = A^L-^(lnL)-P™ , (38) 



for L > 24. The obtained estimate for pm is listed in Table |T|. For the random-bond four 
state Potts model we find Pm = .082(13) while for the models Co_2 the error in pi is at least 
of the same order as the estimate for p^ itself. 

3. Polarization and Xp 

Since the polarization is an order parameter of the transition, the critical polarization, 
P, and Xp should follow the forms 

P ~ L^ , (39a) 

Xp-L- . (39b) 

The critical exponents u, (3p and 7^ are known analytically p^j2^ all along the critical 
line of the pure AT model. For example the ratio — varies continuously between its value 
at the decoupled Ising critical point, ^ = |, and its value at the four state Potts critical 
point, ^ = |- These two extreme values can be obtained, through their relations with ^, 
for any ( pure, anisotropic or disordered) decoupled Ising or four-state Potts model. For the 
decoupled Ising models A^ = A^ = (or A = for the pure model), the a spins system and 
the r spins system are decoupled, hence 

P = (ar) = ((T)(r) = M2~t2/3_ ^4q) 

This implies j3p = 2[3 for all the decoupled Ising models. For the four state Potts models 
K^ = A^ and K^ = A^ (or K = A for the pure four state Potts) so that a spins, r spins and 
ar spins are equivalent, and 

P = (ar) = (a) = (t) = M . (41) 
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It follows that Pp = p. Thus interpreting our results as ^ = | for the critical plane of the 
RBAT implies ^ = 2^ = ^ for the random-bond decoupled Ising model, and ^ = ^ = \ for 
the random-bond four-state Potts model. The question still remains whether the variation 
between these two values along the critical plane of the RBAT will be continuous, as is 
the case for the pure model, or whether it will become discontinuous under the influence of 
disorder. The results presented below are in favor of a continuous variation. The values of 
^ can be obtained through the scaling relation, 

2^ + l^ = d, (42) 

V V 

which our results fully confirm. This leads to ^ = | for the Ising models and ^ = | for the 
four-state Potts models. 

Figures |TT] and |T^ show the critical Xp aiid polarization, P, respectively, as a function 
of logL for the five critical models, C0..4, of the AAT. The solid lines are linear fits to the 
forms (|39|) for L > 24, yielding estimates for critical exponents — and — which are listed in 
Table |. As expected we find that — varies between — = 1.501(1) for the decoupled Ising 
model and — = 1.750(4) ^ - for the four state Potts model. — varies accordingly fulfilling 
the scaling relation ( H2| ). 



Figures |I3] and |TJ show the critical Xp ^^^ polarization, P, respectively, as a function 
of logL for the five critical models, Co. .4, of the RBAT. The same analysis has been carried 
out as for the AAT and the estimated critical exponents — and — are listed in Table HI. 



The results for the random-bond four-state Potts, C4, agree with P = M and Xp = X to a 
high precision and thus need no further explanation. 

There is a clear mismatch between the estimated values of — and — and their expected 
values for the random-bond Ising model Cq. Even the left side of the relation (^OD, P = M^, 
which implied — = 2-, is far from being met, but rather we find P > M^ both for the 
anisotropic and for the random-bond Ising models. This is probably connected to the fact 
that we are measuring, as is usual in a MC calculation, (|M|) and {\P\) which are larger than 
lim/i^o(^) and lim/ip^o(-P) respectively. Both in the random-bond and in the anisotropic 
Ising model the relative variance of P is larger than the relative variance of M. This makes 
(|P|) a worse estimator of \im.hp^o{P) than (|M|) is of \imh^o{M). The discrepancy is a 
finite size effect [Q and should decay as L increases. Indeed if we estimate — from lattices 
of size L > 96 we obtain — = .266(22) which is consistent with the exact value — = 7. An 
estimate of — from lattices of size L > 64 yields — = 1.524(23) which is consistent with the 
exact value — = 1. 

The results for Xp and P seem to favor a continuous variation of — and — for the RBAT 
models as is the case for the AAT models. As explained above, for symmetry reasons one 
expects that if Pismg = /^Potts then /3pjsing 7^ /^p.Potts- This is consistent with the variation of 
exponents we find. None the less since we found the specific heat of all models to diverge 
logarithmically one could expect the RBAT models to be more universal. A possible scenario 
to expect could be Pp = /3p Potts for all RBAT models but the Ising model. Then assuming 
{P/j^)ismg = 1/8 with (H) would imply ^ = 3/2 for the Ising model and ^ = 7/4 for all the 
other RBAT models. Such a behaviour would of course be masqueraded by crossover, but 
none the less we do not find evidence for such a behavior in our results. Further investigation 
of this problem is clearly needed. 
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V. SUMMARY AND DISCUSSION 

We've examined the effect of bond-disorder on tlie line of critical points of the AT model. 
A duality transformation was used in order to locate a critical manifold of the random bond 
AT model which corresponds to the line of critical points of the pure model. An anisotropic 
AT model was used as a convenient reference pure model. Our results consist of the finite 
size dependence of C, M, P, x, Xp ^^ criticality of several random bond and anisotropic AT 
models. The models Co.. .4 interpolate between Ising and four-state Potts models with large 
disorder, thus also monitoring the effect of changing the crossover exponent (a/z/)p. The 
models A2,B2, C2 serve to monitor the effect of changing the amount of disorder while 
keeping the crossover exponent approximately the same ( .37 < {a/u)p < .40 ). 

The results of the critical specific heat of all the RBAT models agree very well with the 
crossover behaviour (^), according to which for L < Lc, C diverges with the same exponent 
(a/z/)p as the corresponding AAT model, while for L > L^, C grows as logL. Our estimates 
for the crossover lengths of the models A2, B2 are much smaller than estimates of Wang et 
al. for random bond Ising models with the same amount of disorder. This is consistent with 
the Harris criterion and the crossover exponent (a/z/)p ~ .4 of the former models and the 
crossover exponent {a/u)^ = for the latter models. For the random bond Ising model we 
find a double logarithmic divergence (PT|), or (^) asymptotically, as was found by Wang 
et al. We view the adequacy of (pS]) for all the RBAT models ( with different (a/z/)p and 
different amounts of randomness) for the full lattice range < L < 256, together with the 
consistency of the crossover lengths with physical understanding (Harris criterion), as a 
strong indication that the RBAT models ( including the four state Potts model ) exhibit a 
a logarithmic divergence of C. We tend to rule out the possibility of a double logarithmic 



divergence of C as predicted by DD |^ for the Baxter model, for reasons explained in sec. 

m 

Although P/u and ■y/u do not vary for the pure model, our results indicate a possibility 
that P/u and 7/z/ do vary for the random models. None the less this very small variation 
might be an artifact of a correction to scaling, e.g. of a logarithmic type. The exponents 
connected with the polarization P and the susceptibility Xp? Pp/^ and 7p/z^ which vary 
continuously for the pure AT model seem to do so also for the RBAT model. But on the 
basis of universality the possibility that these exponents are actually the same (with some 
non-universal corrections to scaling and a crossover effect ) for all RBAT models excluding 
the Ising model, seems to be a plausible scenario in need of further investigation. 

Very recently an experimental study of a two dimensional four-state Potts system with 
quenched impurities was done [^ and the values of z/ = 1.03(8), P = .135(1), and 7 = 
1.68(15) were obtained. These are consistent with our findings for the four-state Potts 
system C4 but with a value of /3 which is smaller than ours. 
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APPENDIX: THE DUALITY TRANSFORMATION OF THE AT MODEL 

We list here the duahty transformation D[K, A] of the AT model and a few of its proper- 
ties. For convenience we use instead of [K, A) different variables for the coupling constants, 

Z = exp-4^ X = exp-2(^+^) X = (X, Z). (Al) 

Under a duality transformation D, a bond of strength X = {X, Z) transforms into a bond 
on the dual lattice of strength X = (X, Z) = D{X) given by [QQ 

X = (1 - Z)/A 

Z= (1-2X + Z)/A. (A2) 

A = 1 + 2X + Z 

We list several properties of the transformation. It is self-inverse, i.e., D[D{X)] = X. It 
maps the zero temperature point ( X = (0, 0) ) onto T = cxd ( X = (1, 1) ) and vice versa. 
The Ising subspace Z = X^ and the 4-state Potts subspace Z = X are invariant under D. 
The subspace S, defined by X > Z > X^ ( or ii" > A > ), is also invariant under D. The 
line Z = 1 — 2X is a self dual line and its intersection with S is the line of critical points of 
the AT model. It is also true that ii Z < 1 — 2X then Z > 1 — 2X and vice versa. 
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FIGURES 

FIG. 1. The original lattice G is denoted by solid lines, while its sites are denoted by dots. The 
dual lattice D is denoted by dashed lines, while its sites are denoted by crosses. A bond of strength 
K is denoted by a thick solid line, while its dual of strength K is denoted by a thick dashed line. 

FIG. 2. Specific heat, C, as a function of logL for five critical models, Co.. 4, of the AAT. Note 
that Co is an anisotropic decoupled Ising model and C4 is an anisotropic four-state Potts model. 
The curves are fits to ( ^) , yielding estimates for ^ which are listed in table | . 

FIG. 3. Specific heat, C, as a function of logL for five critical models, Co. .4, of the RBAT. 
Note that Co is a random-bond decoupled Ising model and C4 is a random-bond four-state Potts 
model. The curves are fits to the form ( EtI) for latice sizes L > 24. 



FIG. 4. Specific heat, C, as a function of logL for five critical models, Co, .4, of the RBAT. 
The curves are fits to the form ( ^8|) , yielding estimates for the coefficients of ( |28D which are listed 
in table . 

FIG. 5. Specific heat, C, as a function of logL for three critical models,742,i?2) C2, of the 
Random Bond AT . The curves are fits to the form ( ^8|) , yielding estimates for the coefficients of 
( ^8|) which are listed in table || . 

FIG. 6. Specific heat, C, as a function of logL for the RBAT critical models, i?2 and Co. ..4 
with fitting curves of ( |3^) . 

FIG. 7. Specific heat, C, as a function of logL for the RBAT critical models, i?2 and Co. ..4 
with fitting curves of ( |3^) . 

FIG. 8. Susceptibility as a function of log L for five critical models, Co. .4, of the AAT. For the 
sake of clarity x ^or the model Cj has been multiplied by 2*. The solid lines are linear fits to the 
form ( ^5|) for L > 24, yielding estimates for critical exponents ^ which are listed in table | . 

FIG. 9. Susceptibility as a function of log L for five critical models, Co.. 4, of the RBAT. For 
the sake of clarity x for the model Cj has been multiplied by 2*. The solid lines are linear fits to 



the form ( ^ ) for L > 24, yielding estimates for critical exponents - which are listed in table III 



FIG. 10. Magnetization as a function of logL for five critical models, Co. .4, of the RBAT. For 
the sake of clarity M for the model Cj has been multiplied by 2*. The solid lines are linear fits to 



the form ( B^ ) fo^ ^ — ^^' yielding estimates for critical exponents - which are listed in table III 
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FIG. 11. Xp ^s a function of logL for five critical models, Co. .4, of the AAT. The solid lines 
are linear fits to the form ( |39b| ) for L > 24, yielding estimates for critical exponents — which are 
listed in table | . 

FIG. 12. Polarization P as a function of logL for five critical models, Co. .4, of the AAT. The 
solid lines are linear fits to the form ( |39a| ) for L > 24, yielding estimates for critical exponents ^ 
which are listed in table | . 

FIG. 13. Xp ^-s a function of logL for five critical models, Co. .4, of the RBAT. The solid lines 
are linear fits to the form ( |39b| ) for L > 24, yielding estimates for critical exponents ^ which are 



listed in table III . 



FIG. 14. Polarization P as a function of logL for five critical models, Co. .4, of the RBAT. The 
solid lines are linear fits to the form ( |39a| ) for L > 24, yielding estimates for critical exponents — 



which are listed in table III . 
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TABLES 

TABLE L Critical exponents ratios and fitting parameters from tlie first series of simulations, 
Co. .4, and the second series A,B,C2 of the Anisotropic AT model. The exponent ratios and 
parameters are defined in the text. 





a 


1 


/3 


> 


Pv 




V 


V 


y 


V 


V 


Co 


.0001(150) 


1.7511(7) 


.1244(4) 


1.501(1) 


.2494(6) 


Ci 


.171(5) 


1.751(1) 


.1242(5) 


1.554(1) 


.2231(5) 


C2 


.375(5) 


1.755(1) 


.1223(7) 


1.608(2) 


.197(1) 


C3 


.549(8) 


1.759(4) 


.1199(27) 


1.667(5) 


.169(3) 


C4 


.630(8) 


1.749(4) 


.1238(21) 


1.750(4) 


.1237(26) 


B2 


.371(5) 


1.753(1) 


.1233(7) 


1.608(2) 


.197(1) 


A2 


.40(1) 


1.750(1) 


.1254(8) 


1.603(2) 


.199(1) 



TABLE n. Fitting parameters for the fits in figures ^ and g according to equations ( |2^ 
Errors are given in parentheses only when error is smaller than or of the same order as the number 
itself. 





ao 


&o 


Co 


Lc 


Co 


-.37(12) 


.58(1) 


5.2E4(1.5E4) 


1. 


Cx 


-4.6 


.51(2) 


1.5E6 


1. 


C2 


-4.1 


.46(127) 


5.5E4 


1. 


C3 


-3.9 


.43(4) 


5.5E4 


1. 


Ca 


-4.1 


.42(1) 


1.0E5 


1. 


B2 


-.09(5) 


2.00(4) 


1.47(10) 


4.0(4) 


A2 


-.07(6) 


9.35(33) 


.26(2) 


51(7) 



TABLE in. Critical exponents ratios and fitting parameters from the first series of simulations, 
Co. .4, and the second series A,B,C2 of the Random Bond AT. These exponents and parameters 
related to the scaling of the order parameters are defined in the text. 





1 

V 


Vx 


V 


■Pm 


> 

V 


V 


Co 


1.751(5) 


-.007(20) 


.125(3) 


-.003(12) 


1.549(9) 


.227(5) 


Ci 


1.751(6) 


-.008(19) 


.124(3) 


-.004(11) 


1.575(8) 


.214(5) 


C2 


1.743(5) 


.026(22) 


.129(3) 


.015(14) 


1.597(9) 


.205(5) 


C3 


1.736(3) 


.057(14) 


.133(2) 


.032(9) 


1.638(5) 


.185(3) 


C4 


1.714(5) 


.148(21) 


.145(3) 


.082(13) 


1.714(5) 


.145(3) 


B2 


1.738(4) 


-.049(17) 


.132(3) 


.029(10) 


1.586(6) 


.209(4) 


A2 


1.739(5) 


-.042(22) 


.132(3) 


.027(13) 


1.590(8) 


.208(5) 
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